program example_MySparseOperator
    use MySparseOperator
    use common_functions
    implicit none
    real(kind=8), dimension(:), allocatable :: a
    integer, dimension(:), allocatable :: rowind, colind
    real(kind=8), dimension(:), allocatable :: b, x

    type(mso_dsparse_matrix) :: Amtx, Bmtx, Cmtx

    integer :: msglvl
    integer :: maxiter(2), itercount
    real(kind=8) :: maxrelres, residual(2)
    integer :: nnz

! preallocate the sparse matrixs
    nnz = 40
    call mso_sparse_d_preallocate(Amtx,9,9,nnz)
    write(*,*) size(Amtx%ia),size(Amtx%ja),size(Amtx%a)
    call mso_sparse_d_preallocate(Bmtx,9,9,nnz)
    call mso_sparse_d_preallocate(Cmtx,9,9,nnz)

! create a sparse matrix by indices of elements
    write(*,*) '============= create sparse matrixs by indices of elements =============='
    call mso_sparse_d_create_coo( (/8:1:-1/), (/ 8:1:-1 /), (/ real((/1:8/),8)/), Amtx, 8,10)
    call mso_sparse_show(Amtx, 'f5.1')
    call mso_sparse_d_create_coo( (/8:1:-1/), (/ 8:1:-1 /), (/ real((/1:8/),8)/), Amtx, 10,8)
    call mso_sparse_show(Amtx, 'f5.1')
    write(*,*) size(Amtx%ia),size(Amtx%ja),size(Amtx%a)

    call mso_sparse_d_create_coo( (/9:1:-1/), (/ 9:1:-1 /), (/ real((/1:9/),8)/), Amtx, 9)
    call mso_sparse_show(Amtx, 'f5.1')
    write(*,*) size(Amtx%ia),size(Amtx%ja),size(Amtx%a)

! add two sparse matrix
    write(*,*) '============= add two sparse matrix ======================================'
    call cf_allocate(a,9)
    a = 3.d0
    call mso_sparse_d_create_coo( (/1:8,1/), (/2:9,1/), a, Bmtx, 9)
    call mso_sparse_d_add(Amtx,1.d0,Bmtx,Cmtx)
    call mso_sparse_show(Cmtx, 'f5.1')
    call mso_sparse_show(Cmtx, 'f5.1','n',-1)
    write(*,*) size(Cmtx%ia),size(Cmtx%ja),size(Cmtx%a)
   write(*,*) '=========================================================================='
    call mso_sparse_d_add(Cmtx,1.d0,Cmtx,Amtx,'t')
    call mso_sparse_show(Amtx, 'f5.1')

! delete two rows and two columns from a sparse matrix
    write(*,*) '====== delete two rows and two columns from a sparse matrix =============='
    write(*,*) '=== delete two rows =========================='
    write(*,*) '=== before deleting ================='
    call mso_sparse_show(Amtx,'f5.1')
    write(*,*) size(Amtx%ia),size(Amtx%ja),size(Amtx%a)
    call mso_sparse_d_deleteRowCol(Amtx,(/2,3/),'R',Bmtx)
    write(*,*) '=== after deleting ================='
    call mso_sparse_show(Bmtx, 'f5.1')

    write(*,*) '=== delete two columns ======================'
    write(*,*) '=== before deleting ================='
    call mso_sparse_show(Bmtx,'f5.1')
    write(*,*) size(Bmtx%ia),size(Bmtx%ja),size(Bmtx%a)
    call mso_sparse_d_deleteRowCol(Bmtx,(/5,7, 10, 11/),'C',Amtx)
    write(*,*) '=== after deleting ================='
    call mso_sparse_show(Amtx, 'f5.1')

! a sparse matrix multiply a vector
    write(*,*) '====== a sparse matrix multiply a vector =============='
    call mso_sparse_create_coo( (/1:9,1:8,2:9/), (/1:9,2:9,1:8/), real((/1:9,1:8,2:9/),8), Amtx, 9)
    call cf_allocate(x,9)
    call cf_allocate(b,9)
    x = 1.d0
    call mso_sparse_d_mv(Amtx,x,b)
    write(*,*) '===== the matrix is ============'
    call mso_sparse_show(Amtx,'f5.1')
    write(*,*) '===== the vector is ============'
    call prtVec(x,'f5.1')
    write(*,*) '===== the result is ============'
    call prtVec(b,'f5.1')

! two sparse matrix multiply with each other
    write(*,*) '====== two sparse matrix multiply with each other =============='
    call mso_sparse_create_coo( (/1:9,1:8,2:9/), (/1:9,2:9,1:8/), real((/1:9,1:8,2:9/),8), Amtx, 9,9)
    call mso_sparse_create_coo( (/2,3,5,6,7,8/), (/1:6/), real((/1,1,1,1,1,1/),8), Bmtx, 9,6)
    write(*,*) '====== the first sparse matrix A is ======'
    call mso_sparse_show(Amtx,'f5.1')
    write(*,*) '====== the second sparse matrix B is ======'
    call mso_sparse_show(Bmtx,'f5.1')
    write(*,*) '====== C=A*B ======'
    call mso_sparse_mm(Amtx,Bmtx,Cmtx)
    call mso_sparse_show(Cmtx,'f5.1')
    write(*,*) size(Cmtx%ia), size(Cmtx%ja), size(Cmtx%a)
    write(*,*) '====== B^T*C= ======'
    call mso_sparse_mm(Bmtx,Cmtx,Amtx,'t')
    call mso_sparse_show(Amtx,'f5.1')
    write(*,*) size(Amtx%ia), size(Amtx%ja), size(Amtx%a)



! get the upper triangular part of a sparse matrix
    write(*,*) '========= get the upper triangular part of a sparse matrix ==============='
    write(*,*) '====== upper =  0 ==========='
    call mso_sparse_uppertriangle(0,Amtx,Bmtx)
    call mso_sparse_show(Bmtx, 'f5.1')
    write(*,*) '====== upper =  1 ==========='
    call mso_sparse_uppertriangle(1,Amtx,Bmtx)
    call mso_sparse_show(Bmtx, 'f5.1')
    write(*,*) '====== upper =  -1 ==========='
    call mso_sparse_uppertriangle(-1,Amtx,Bmtx)
    call mso_sparse_show(Bmtx, 'f5.1')


! use LU method (mso_dpardiso) to caculate the solution of a sparse linear system
    write(*,*) '======= caculate the solution of a sparse linear system ================='
    call mso_sparse_d_create_coo( (/ 1:9 /), (/ 1:9 /), 10.d0*(/1:9/), Amtx, 9)
    call mso_sparse_d_create_coo( (/ 1:8/), (/ 2:9/), 1.d0*(/1:8/), Bmtx, 9)
    call mso_sparse_add(Amtx,1.d0,Bmtx,Cmtx)
    call mso_sparse_add(Cmtx,1.d0,Cmtx,Amtx,'T')
    call mso_sparse_copy(Amtx,Cmtx)
    call cf_allocate(x,9)
    x = real((/1:9/),8)
    call mso_sparse_mv(Cmtx,x,b)
    write(*,*) '======= the matrix is:'
    call mso_sparse_show(Cmtx,'f7.1')

    write(*,*) '======= expected solution is:'
    call prtVec(x,'f5.1')
    call mso_sparse_d_pardiso(Cmtx,b,x)
    write(*,*) '======= computed solution is:'
    call prtVec(x,'f5.1')


! another example to use LU method (mso_dpardiso) to solve sparse linear system
    write(*,*) '======== another example about using LU method (mso_dpardiso) to solve sparse linear system ======'
    call cf_allocate(rowind,18)
    rowind = (/ 1, 1,          1, 1,    &
                   2, 2,    2,          &
                      3,             3, &
                         4,       4,    &
                            5, 5, 5,    &
                               6,    6, &
                                  7,    &
                                     8  /)
    call cf_allocate(colind,18)
    colind = (/ 1, 2,          6, 7,    &
                   2, 3,    5,          &
                      3,             8, &
                         4,       7,    &
                            5, 6, 7,    &
                               6,    8, &
                                  7,    &
                                     8  /)
    call cf_allocate(a,18)
    a = (/ 7.d0,  1.d0,                    2.d0,  7.d0,        &
                 -4.d0, 8.d0,       2.d0,                      &
                        1.d0,                           5.d0,  &
                              7.d0,               9.d0,        &
                                    5.d0,  1.d0,  5.d0,        &
                                          -1.d0,        5.d0,  &
                                                 11.d0,        &
                                                        5.d0   /)
    call mso_sparse_create_coo(rowind,colind,a,Amtx,8)
    call mso_sparse_add(Amtx,1.d0,Amtx,Bmtx,'T')
    call cf_allocate(x,8)
    x = 1.d0    ! expected solution
    call cf_allocate(b,8)
    call mso_sparse_mv(Bmtx,x,b)
    print *, 'the matrix is:'
    call mso_sparse_show(Bmtx,'f6.1')
    print *, 'the right hand:'
    call prtVec(b,'f5.1')
    print *, 'expected solution:'
    call prtVec(x,'f5.1')
    call mso_sparse_setype(Bmtx,-2)
    call mso_sparse_pardiso(Bmtx,b,x)
    print *, 'computed solution:'
    call prtVec(x,'f5.1')


! use GMRES to solve example above
    write(*,*) '======== use GMRES to solve sparse linear system ========================='
    maxiter(1) = 20 ! number of outer iteration
    maxiter(2) = 9 ! number of inner iteration
!     read(*,*) maxiter
    maxrelres = 1.d-8 ! specifies the relative tolerance
    x = 1.d0
    print *, 'expected solution:'
    call prtVec(x,'f5.1')
    call mso_sparse_mv(Bmtx,x,b)
    x = 0.d0 ! initial solution
    msglvl = 1 ! optional, print the data of computation; else, set to be 0, not shown any information
    call mso_sparse_d_gmres(maxiter, maxrelres, Bmtx, b, x, itercount, residual,msglvl)
    write(*,*) 'itercount=', itercount
    write(*,*) 'absolue residual=', residual(1)
    write(*,*) 'relative residual=', residual(2)
    print *, 'computed solution:'
    call prtVec(x,'f5.1')




    call cf_deallocate(a)
    call cf_deallocate(b)
    call cf_deallocate(x)
    call cf_deallocate(rowind)
    call cf_deallocate(colind)
    call mso_sparse_destroy(Amtx)
    call mso_sparse_destroy(Bmtx)
    call mso_sparse_destroy(Cmtx)

end program example_MySparseOperator
